function [gauss_points, gauss_weights] = gauss_points_weights_tri(N_gauss_int2d)

switch N_gauss_int2d
    case 3
        gauss_points = [0.5, 0;
            0, 0.5;
            0.5, 0.5];
        gauss_weights = [1/3, 1/3, 1/3];
    case 4
        gauss_points = [1/3, 1/3;
            0.6, 0.2;
            0.2, 0.6;
            0.2, 0.2];
        gauss_weights = [-27/48, 25/48, 25/48, 25/48];
    case 6
        gauss_points = [0.816847572980459, 0.091576213509771;
            0.091576213509771, 0.816847572980459;
            0.091576213509771, 0.091576213509771;
            0.108103018168070, 0.445948490915965;
            0.445948490915965, 0.108103018168070;
            0.445948490915965, 0.445948490915965];
        gauss_weights = [0.109951743655322, 0.109951743655322, 0.109951743655322, 0.223381589678011, 0.223381589678011, 0.223381589678011];
    case 9
        gauss_points = [1/2, 1/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))*(1+sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))*(1-sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))*(1+sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))*(1-sqrt(3/5))/4;
            1/2, (1+sqrt(3/5))/4;
            1/2, (1-sqrt(3/5))/4;
            (1+sqrt(3/5))/2, (1-sqrt(3/5))/4;
            (1-sqrt(3/5))/2, (1+sqrt(3/5))/4];
        gauss_weights = [16/81, 200/324*(1-sqrt(3/5))/8, 200/324*(1-sqrt(3/5))/8, 200/324*(1+sqrt(3/5))/8, 200/324*(1+sqrt(3/5))/8, 10/81, 10/81, 10/81*(1-sqrt(3/5)), 10/81*(1+sqrt(3/5))];
    case 12
        gauss_points = [0.873821971016996, 0.063089014491502;
            0.063089014491502, 0.873821971016996;
            0.063089014491502, 0.063089014491502;
            0.501426509658179, 0.249286745170910;
            0.249286745170910, 0.501426509658179;
            0.501426509658179, 0.501426509658179;
            0.636502499121399, 0.310352451033785;
            0.636502499121399, 0.053145049844816;
            0.310352451033785, 0.636502499121399;
            0.310352451033785, 0.053145049844816;
            0.053145049844816, 0.636502499121399;
            0.053145049844816, 0.310352451033785];
        gauss_weights = [0.050844906370207, 0.050844906370207, 0.050844906370207, 0.116786275726379, 0.116786275726379, 0.116786275726379, 0.082851075618374, 0.082851075618374, 0.082851075618374, 0.082851075618374, 0.082851075618374, 0.082851075618374];
    otherwise
        error('Number of quadrature point should be 3,4,6,12.')
end

end